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Abstract 

Marching  techniques  for  the  parabolized  Navier  Stokes  equa¬ 
tions  are  considered.  With  the  full  pressure  interaction  and 
prescribed  edge  pressure  these  equations  are  weakly  elliptic  in 
subsonic  zones.  A  minimum^  marching  step  size  (Axm^n) >  Pro¬ 
portional  to  the  total  thickness  (y,.)  of  the  subsonic  layer, 
exists.  However,  for  thin  subsonic  boundary  layers  (yM  <<  1) 
and  with  Ax  =  0 !( yM)  ,  stable  and  accurate  solutions  are  possible. 
With  forward  differencing  of  the  axial  pressure  gradient  the 
procedure  can  be  made  unconditionally  stable?  a  global  iteration 
procedure,  requiring  only  the  storage  of  the  pressure  term,  has 
been  demonstrated  for  a  separated  flow  problem.  Solutions  for 
incompressible  boundary  layer-like  flows,  for  internal  flows, 
and  for  supersonic  flow  over  a  cone  at  incidence  with  a  coupled 
strongly  implicit  procedure  are  presented. 

1.  Introduction 

It  has  now  been  generally  accepted  that  boundary  layer 
methodology  can  be  extended  to  the  so-called  parabolized  Navier- 
Stokes  (PNS)  equations  for  a  significant  variety  of  flow  problems. 
In  a  recent  paper,  Davis  and  Rubin^  have  reviewed  several  viscous 
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flow  computations  in  which  parabolized  or  thin  layer  techniques 
have  been  applied  in  order  to  accurately  determine  the  flow 
characteristics.  This  publication  also  reviews  some  of  the 
early  history  of  the  PNS  development. 

The  purpose  of  the  present  paper  is  to  discuss  some  recent 
investigations  using  the  PNS  equations.  In  particular,  we  are 
concerned  here  with  efficient  three-dimensional  algorithms,  a 
clearer  understanding  of  the  limits  of  applicability  of  PNS 
marching  techniques,  and  pressure  interaction  relaxation  for 
separated  flows  and  other  problems  where  upstream  influence  is 
of  importance.  In  this  regard,  three  solution  procedures  are 
considered:  (1)  single-sweep  "boundary  layer-like"  marching  for 

two  and  three-dimensional  flows;  (2)  multiple  sweep  iteration  or 
global  pressure  relaxation,  where  upstream  influence  and  possibly 
axial  flow  separation  are  important,  but  regions  of  subsonic 
flow  are  small;  and  (3)  global  relaxation  where  subsonic  flow 
domains  are  large.  For  the  latter  two  classes  of  problems,  the 
analysis  draws  heavily  on  that  of  interacting  boundary  layers 
and  inviscid  subsonic  relaxation  methods  where  applicable. 

In  the  course  of  proceeding  to  specific  examples,  a  brief 
review  of  the  limitations  associated  with  PNS  marching  or  relaxation 
is  necessary.  The  PNS  equations,  which  for  simplicity  are  given 
here  only  for  "two-dimensional  incompressible  flow,"  are  as 
follows: 


u.  +  uu  +  vu  = 
t  x  y 


Px  +  R  Uyy  ' 


(la) 


2 


(lb) 


V,  +  UV  +  VV  =  -  P  +  (i  V  )  , 

t  x  y  *y  R  yy 

Ux  +  Vy  =  0-  (lC) 

R  is  the  Reynolds  number;  x  is  denoted  the  axial  flow  direction 

and  y  the  normal  direction.  The  system  (1)  differs  from  the 

complete  Navier-Stokes  equations  only  by  the  omission  of  the 

axial  diffusion  terms.  Strictly  speaking  the  inclusion  of  the 

v  term  in  (lb)  is  inconsistent  with  the  omission  of  u 
yy  xx 

in  (la) .  Either  the  former  should  be  neglected  (this  is  probably 
a  more  appropriate  definition  of  the  PNS  equations)  or  the 
latter  retained.  In  fact,  these  terms  have  little  effect  on 
any  of  the  results  presented  here.  The  mathematical  character 
of  (1)  is  controlled  by  the  px  term  in  (la) .  When  px  is  pre¬ 
scribed  (assumed  known) ,  the  system  is  parabolic.  This  was  the 

2 

case  in  the  original  merged  layer  analysis  of  Rudman  and  Rubin  , 
where  for  hypersonic  cold  wall  flow  px  can,  in  fact,  be  neglected 
in  (la).  It  should  be  emphasized,  however,  that  axial  pressure 
gradients  are  still  present  and  are  evaluated  through  the  momentum 
equation  (lb) ,  and  the  energy  equation  in  the  compressible  case. 

When  the  px  term  in  (la)  is  retained  implicitly  (not 
prescribed) ,  the  system  (1)  is  no  longer  parabolic,  as  an  elliptic 
pressure  or  acoustic  interaction  occurs  in  regions  of  subsonic 
flow.1  The  "parabolic"  form  for  the  velocities  has  led  to  the 
expression  parabolized  Navier-Stokes  equations.  This  pressure 
interaction  also  appears  for  boundary  layer  equations,  when  p 
is  not  prescribed.  The  resulting  upstream  interaction  has  been 
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analyzed  by  Lighthill5,  who  demonstrated  the  existence  of 
exponential  growing  solutions  in  the  interaction  zone.  Similar 
behavior  has  been  encountered  with  marching  procedures  for  the 
PNS  equations.  The  primary  difference  between  the  boundary 
layer  and  PNS  equations  is  that  for  the  latter  the  pressure 
interaction  is  manifested  through  both  the  outer  pressure 
boundary  formulation  and  the  normal  momentum  equation  (lb) . 

2.  Single  Sweep  Marching 

For  problems  where  upstream  influence  and  axial  flow  separation 
are  not  significant,  it  is  natural  to  consider  the  system  (1) 
by  boundary  layer  marching  techniques,  i.e.,  backward  differences 
are  applied  for  all  x-gradients.  If  p  is  prescribed,  this 

A 

approach  is  quite  acceptable  as  the  equations  are  in  fact 
parabolic.  For  implicit  numerical  schemes,  the  marching  calcu¬ 
lation  should  be  unconditionally  stable  for  all  Ax  marching 
steps,  see  Davis  and  Rubin1  for  additional  references.  On  the 
other  hand,  if  px  is  assumed  unknown,  the  "elliptic"  pressure 
interaction  of  Lighthill  is  introduced  and  therefore  the  exponential 

growth,  representative  of  upstream  influence,  can  be  anticipated. 

4 

Lubard  and  Helliwell  have  examined  the  stability  of  the  backward 

difference  approximation  for  px  in  (la)  and  they  have  shown  that 

for  Ax  <  (Ax)m^n  instability  or  departure  solutions  will  occur. 

Similar  results  were  found  earlier  by  Lin  and  Rubin.5  For 

Ax  >  (Ax)  .  the  marching  scheme  is  stable.  Therefore  (Ax)  . 

mm  mm 

would  appear  to  represent  a  measure  of  the  upstream  elliptic 
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interaction.  For  Ax  <  (Ax)  .  the  marching  scheme  attempts  tc 

mm 

represent  this  interaction  and  therefore  the  Lighthill  behavior 
should  be  recovered. 

Since  the  backward  difference  formula  for  px  does  not  provide 

any  upstream  contribution,  it  does  not  properly  represent  the 

differential  form  for  Ax  <  (Ax)  .  .  When  the  backward  difference 

mm 

approximation  is  less  representative  of  px,  the  error  introduced 

serves  to  reduce  (Ax)  For  example,  at  x  =  x.,  (p.  n  -p.  ~)/Ax 

is  less  severe  than  (p^  -  p^_^) /Ax,  and  for  px»  prescribed, 

(Ax)  .  =0.  In  view  of  this  behavior,  several  investigators, 

mm  '  ^ 

6  7  8 

Vigneron  et  al.  ,  Yanenko  et  al.  and  Lin  and  Rubin  have  attempted 

to  eliminate  the  pressure  interaction  by  incorporating  "small" 

inconsistencies  into  the  difference  approximation.  They  have 

assumed  (1)  a  variable  Ax  in  the  difference  form  for  px  such 

that  Ax  >  (Ax)min  locally6,  (2)  "regularization"  functions  of 

the  type  (1//R)  f(ux'vx»Px)  as  modified  coefficients  for  the 

7 

uux  and  px  terms  in  (la)  ,  and  (3)  the  use  of  finite  temporal 
iteration  to  modify  the  convection  velocity  in  each  step  of 

g 

the  marching  procedure.  Each  of  these  techniques  introduces 
some  inconsistency  into  the  difference  equations  in  subsonic 
zones;  reasonable  results  have  been  obtained  with  these  methods 
for  certain  problems.  In  order  for  these  techniques  to  be 
effective,  the  inconsistency  must  be  large  enough  to  suppress 
the  elliptic  character,  yet  small  enough  to  maintain  an  acceptable 
order  of  accuracy. 

The  PNS  model  has  been  considered  in  some  detail  by  Rubin 

g 

and  Lin.  For  the  system  (1)  with  a  backward  difference  for 
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3/9x  and  central  y  differences,  the  linear  von  Neumann  stability 
analysis  leads  to  the  following  condition  for  the  eigenvalues  X: 


2  8 

u(X-l)  + 4b  Xsin  ^  +  I  c  XsinB 
u(X-l)  cos^  |-  +  4b  Xsin2  j  +  I  c  Xsing 

2 

where  a  =  Ax/Ay,  b  =  Ax/RAy  ,  c  =  av,  I  =  /-l,  F  =  ( X— 1) . 

It  can  be  shown  that  the  value  of  X  is  closely  related 

max 

to  the  highest  frequency  mode,  so  that  when  the  number  of  grid 
points  across  the  layer  N  >>  1,  3  =  tt/(N-1).  Equation  (2) 
then  takes  the  simplified  form 


(  X-l)  F  ( X) 


4a  X 


.  2  6 
sin  2 


(2) 


(X-l) 


2  2 


M=1, 


(3) 


where  A  =  trAx/yM  and  yM  is  the  layer  thickness?  yM  =  (N)  Ay. 

The  condition  (3)  indicates  that  Ax/yM  is  the  relevant  stability 
parameter.  From  (3) ,  the  marching  procedure  will  be  stable  for 

A  =  TrAx/yM  >  2 

Therefore 


(Ax) 


min 


f  M 


(4) 


The  complete  numerical  solution  of  (2),  for  all  S,  has  been 

2 

obtained,  and  the  analytic  result  min  i  ^  YM  confirmed. 

The  extent  of  the  elliptic  numerical  interaction  is  of  the  order 

of  the  thickness  of  the  total  layer.  If  the  system  (1)  is  used 

-1/2 

to  solve  boundary  layer  problems,  then  yM  =  0(R  )  and 


-1/2 

therefore  (Ax)min  =  0(R  ).  For  interaction  regions  where 

11  2  —3/8 

triple  deck  structure  is  applicable,  (Ax)  .  =  —  y..  *  0(R  '  ) 

or  the  extent  of  the  Lighthill^  upstream  influence.  This  would 

tend  to  confirm  the  idea  of  a  limited  elliptic  zone  contained 

in  the  PNS  formulation.  For  Ax  >  (Ax)  .  ,  this  elliptic  effect 

mm 

is  suppressed.  When  the  upstream  influence  is  negligible,  this 
inconsistency  should  have  little  effect  on  the  solution.  For 
truly  interactive  flows  the  ellipticity  must  be  retained  and  the 
global  forward-difference  concept  discussed  in  the  next  section 
is  required. 

The  results  for  the  PNS  and  other  "transonic"  equations 

clearly  indicate  that  step  sizes  of  the  order  of  the  subsonic 

—1/2  —3/8 

region,  which  for  supersonic  mainstreams  is  0(R  '  )  or  0(Re  '  ) 

in  a  triple  deck  region,  will  provide  stable  and  accurate 
solutions  for  flows  in  which  upstream  effects  are  not  dominant. 

In  a  later  section,  where  a  strongly  implicit  algorithm  is  intro¬ 
duced  to  obtain  marching  solutions  for  the  supersonic  flow  over  a 

cone  at  incidence,  this  (Ax)  .  dependency  on  yw  will  be  shown 

mm  m 

for  the  compressible  PNS  equations. 

3.  Multiple  Sweep  Marching-Global  Iteration 

If  consistency,  for  Ax  -*■  0,  of  the  difference  formulation  is 
to  be  achieved,  or  if  upstream  influence  is  important  and/or 
separation  occurs,  then  backward  differencing  of  the  px  term 
should  be  rejected.  With  any  form  of  forward  or  central  dif¬ 
ferencing  for  p^,  relaxation  (multiple  marching  sweeps  or  global 
iteration)  is  required.  Three  possibilities  will  be  discussed: 
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(1)  forward  differencing,  (2)  central  differencing,  and  (3)  use 
of  the  Poisson  equation  for  pressure. 

Forward  differencing  for  p^  introduces  the  upstream  value 
pi+l'  Prescribed  in  each  marching  sweep,  and  also  has 

the  advantage  of  including  the  local  pressure  value  p^.  This 
provides  for  coupling  of  the  pressure-velocity  system  (1)  and 
allows  for  a  free  surface  pressure  interaction.  This  is  important 
for  problems  where  axial  flow  separation  occurs. 

At  this  point,  a  few  remarks  concerning  the  role  of  the  px 
term  for  separated  flow  are  relevant.  In  recent  years  there 
has  been  considerable  analysis  of  free  pressure-boundary  layer 
interactions  for  separated  flows.10  There  is  general  agreement 
that,  for  limited  separation  bubbles,  boundary  layer  equations 
can  be  used  to  calculate  such  flows.  Moreover,  the  singularity 
at  separation  does  not  appear  if  the  p* (x)  term  is  not  prescribed 
but  allowed  to  develop.  Inverse  methods,  in  which  the  displacement 
thickness  or  shear  stress  is  prescribed  during  the  marching  step 
and  then  updated  in  subsequent  relaxation  sweeps  have  been 
considered,  as  have  procedures  in  which  temporal  terms  are 
retained  in  order  to  introduce  the  pressure  or  displacement 
thickness  implicitly.  In  these  procedures,  the  p1 (x)  term  is 
replaced  with  an  interaction  expression  (displacement  slope  for 
supersonic  flow10  or  Cauchy  integral  for  subsonic  flow11) .  The 
local  pressure  or  displacement  thickness  then  appears  implicitly 
and  is  coupled  with  the  velocity  evaluation.  In  all  of  these 
interaction  analyses,  those  components  of  the  px  approximation 
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that  introduce  upstream  terms  are  updated  during  the  relaxation 
sweeps . 

From  interacting  boundary  layer  analysis  we  can  then 
conclude  that  if  the  elliptic  character  is  to  be  modeled  consis¬ 
tently,  the  p^  representation  for  the  PNS  equations  should 
introduce  downstream  contributions.  If  the  separation  singularity 
is  to  be  circumvented  in  any  relaxation  sweep,  a  free  pressure 
interaction  through  the  outer  boundary  condition  or  through 
the  y-momentum  equation  (lb)  must  be  introduced.  Forward 
differencing  would  appear  to  satisfy  both  of  these  constraints. 

The  stability  analysis  for  equations  (1)  has  been  extended  in 
reference  9  for  a  variety  of  px  approximations.  For  forward 
differencing  of  px,  the  stability  condition  (4)  is  modified 
solely  by  the  factor  F(X),  such  that  F(X)  =  -X.  From  (4)  it 
can  be  inferred  that  for  8  =  ir/N-1,  this  procedure  is  uncon¬ 
ditionally  stable.  The  stability  curves  for  general  8  values  are 
given  in  reference  9.  Forward  differencing  is  unconditionally 
stable  for  all  8  at  all  R.  It  is  significant,  however,  that  as 
the  convective  velocity  v  increases  both  eigenvalues  asymptote 
to  one.  So  that  stability  is  marginal  when  the  subsonic  region 
is  large. 

During  the  first  sweep  of  the  global  iteration  procedure 
the  value  must  be  prescribed  by  an  initial  guess.  p^+^  can 

be  chosen  as  a  constant  equal  to  the  boundary  condition  at  the 
outer  boundary,  or  the  surface,  or  some  combination  thereof. 

Since  the  variation  of  p  across  the  subsonic  layer  is  small,  any 
of  these  values  will  generally  suffice. 


In  order  to  test  the  applicability  of  forward  differencing 
for  p  two  boundary  layer  problems  have  been  considered  with  the 
full  PNS  system  (1) : 


(i)  Flat  Plate:  u=p=l  at  y=  y„; 

M 

u  =  v  =  0  at  y  =  0. 


(ii)  Separation  Bubble:  u=p=l,  x  <  0 


u  =  1-x,  p  =uu,0<x<0.25 

X  X  —  — 


u  =  0.75,  p^  =  0,  x  >_0.25 


} 


at 

y = 


u  =  v  =  0  at  y  =  0 . 


The  pressure  gradient  p  is  forward  differenced.  All  other  x 
derivatives  are  backward  differenced.  These  are  neglected  in 
regions  of  reverse  flow.  All  y  derivatives  are  central 
differenced,  except  for  the  continuity  equation  (lc)  and  momentum 
equation  (lb)  where  the  trapezoidal  rule  is  used.  Multiple 
sweeps  or  global  iteration  was  stable  and  converged  for 
(Ax)  .  2  (y„/6)  ;  the  value  of  (Ax)  .  =  y,./60  was  also  tested 
and  with  forward  differencing  was  stable.  It  is  significant  that 
in  this  relaxation  procedure  it  is  necessary  to  store  only  the 
pressure  field  for  each  successive  iteration  level.  The  velocities 
are  re-evaluated  during  each  marching  sweep.  The  results  are  in 
excellent  agreement  with  published  results  for  both  problems. 

The  free  surface  pressure  interaction  introduced  by  the  y-moemntum 
equation  has  eliminated  the  separation  singularity  for  the  bubble 


problem.  The  value  of  Ax  equal  to  one- sixth  the  boundary  layer 

-1/2 

thickness  yM  =  0(R  )  appears  to  be  adequate.  For  the  smaller 

value  of  (Ax)  .  =  yM/60  or  with  y„  =  (MR-8^8),  i.e.,  the 

min  M  M 

-1/2 

triple  deck  interaction  length,  rather  than  R  the  boundary 
layer  length  scale,  some  variations  in  the  solution  were  obtained. 
These  are  not  given  here.  From  the  stability  results,  we  note 
that  the  value  of  Ax  can  be  made  arbitrarily  small;  however, 
for  the  present  examples,  this  is  unnecessary.  For  the  cone 
geometry,  to  be  considered  in  a  following  section,  considerably 
smaller  values  of  Ax  are  used.  Some  typical  results  for  the 
boundary  layer  examples  are  shown  in  figures  (1)  to  (4) . 


Fig.  1 


Fig.  2 


Fig.  3 


Fig.  4 
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For  the  flat  plate  case,  comparisons  between  the  PNS  and 

3  7 

Blasius  solution  are  shown,  for  R  =  10  and  10  ,  in  figure  (1) 

for  the  velocity  components,  and  in  figure  (2)  for  the  surface 

skin  friction  coefficient  Cf.  The  agreement  is  quite  reasonable. 

The  maximum  error  occurs  at  the  surface  and  this  can  be  seen 

from  the  figures.  Additional  results  for  the  pressure  variation 

across  the  boundary  layer  are  given  in  reference  9.  The  pressure 

p^+1  is  updated  during  each  sweep  of  the  global  iteration 

procedure.  For  the  initial  iteration  p^+^(x,y)  was  taken  equal 

to  the  prescribed  edge  value;  i.e. ,  p(xi+1,y)  =  p(xi+1,yM) . 

During  the  relaxation  process  small  pressure  variations  are 

calculated  across  the  boundary  layer.  The  qualitative  agreement 

g 

with  the  third-order  Blasius  pressure  distribution  is  good. 
Solutions  for  the  separation  bubble  case  are  given  in  figure  (3) 
for  typical  isovels,  and  in  figure  (4)  for  the  surface  pressure 
variation.  The  predicted  separation  point  value  of  xgep  =  0.1180 
is  close  to  the  boundary  layer  value  of  0.1198.  Both  separation 
and  reattachment  points  exhibited  smooth  transitions  and  con¬ 
vergence.  Since  the  outer  boundary  conditions  were  fixed  and  the 
second-order  Cauchy  integral  displacement  condition  was  not 
imposed,  the  free  interaction  was  manifested  solely  through  the 
y-momentum  equation  (lb) .  Inclusion  of  the  displacement  boundary 
condition  should  have  a  slight  effect  on  the  solutions.  Con¬ 
vergence  of  the  global  relaxation  procedure  is  quite  rapid. 

9 

Only  five  to  ten  iterations  are  required.  Of  course,  for  the 
problems  considered  here  the  pressure  variation  across  the  layer 
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is  small  so  that  the  initial  guess  is  quite  good.  In  view  of 
the  stability  analysis  previously  discussed,  and  since  the  PNS 
system  includes  all  of  the  elements  of  both  boundary  layer  and 
triple  deck  equations,  the  present  solutions  with  y^  =  0(R  /  ) 
should  reproduce  the  results  obtained  with  these  approximations. 
Detailed  comparisons  will  be  the  subject  of  future  studies. 

If  central  differencing  is  used  for  p  the  downstream  point 

X 

p^+^  is  introduced  once  again;  however,  the  value  pi  no  longer 
appears  and  the  y-momentum  equation  will  be  uncoupled  from  the 
velocities  unless  p^  is  re-introduced.  Several  possibilities 
exist:  (1)  p^  appears  directly  through  the  outer  boundary 

condition,  as  in  interacting  boundary  layer  theory;10,11 

(2)  a  temporal  relaxation  term  p^.  is  introduced  in  (la)  .  This 

has  been  used  in  ADI  solutions  for  interacting  boundary  layers;10 
12 

(3)  a  K-R  approximation,  where  forward  differencing  is  corrected 
during  the  relaxation  sweeps,  is  applied. 

From  the  stability  analysis  for  each  sweep  of  the  marching 

procedure,  it  is  seen  that  with  central  differencing,  the  function 

F ( X)  =  -1  in  (2).  This  is  an  unconditionally  unstable  condition 

and  therefore  further  reinforces  the  need  for  a  p^  contribution. 

With  an  appropriate  p^  contribution  F(X)  =  (aX-1)  ,  where  a 

reflects  the  p^  term.  From  (2) ,  the  marching  procedure  is 

stabilized  conditionally  for  all  a  ^  0  (recall  the  earlier 

(Ax)  .  condition  for  a  -  1) ;  however,  for  a  <  -1,  unconditional 
min  — 

stability  results.  Actual  experience  with  the  various  px 
approximations  for  the  compressible  PNS  system  is  discussed  in  a 
following  section  on  flow  over  a  cone  at  incidence. 
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Finally,  a  stability  analysis  for  convergence  of  the  global 
iteration  procedure  has  also  been  considered.  When  px  is 
treated  implicitly,  to  some  degree,  preliminary  conclusions  are 
that  convergence  is  assured.  On  the  other  hand,  if  p  is 
treated  explicitly,  i.e.,  from  a  previous  marching  solution, 
it  would  appear  that  the  global  iteration  procedure  will 
diverge.  The  pressure  results  solely  from  the  uncoupled  normal 
momentum  equation  (lb). 

4.  Global  Iteration  for  Subsonic  Flow 

For  problems  where  subsonic  regions  are  not  confined  to  thin 
layers,  single  sweep  backward  differencing  for  px  will  generally 
not  be  reliable  since  the  stability  condition  requires  that 
Ax  >  (Ax)  .  =  0(y„)  =  0(1).  Therefore,  global  relaxation  for 

the  pressure  interaction  is  required  in  all  cases. 

Two  procedures  for  calculating  such  flows  have  been  considered. 

In  the  first,  the  pressure  interaction  is  evaluated  with  the 

Poisson  form  of  the  pressure  equation.  In  the  second,  which  is 

currently  under  development,  the  pressure  is  coupled  directly 

to  the  elliptic  velocity  solver  by  a  splitting  procedure  in  the 

13 

spirit  of  Dodge.  There  are  a  number  of  significant  differences 
however  that  allow  for  a  completely  consistent  global  relaxation 
formulation.  For  both  of  these  formulations  a  Poisson  operator 
appears  explicitly  in  the  equations  and  therefore  the  elliptic 
character  of  the  equations  is  further  strengthened.  Only  the 
former  approach  shall  be  described  here  in  any  detail.  Solutions 
are  presented  for  an  asymmetric  channel  having  a  moderate 
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constriction.  The  complete  analysis  is  given  in  reference  14. 
A  short  review  is  presented  here. 

Equations  (la)  and  (lb)  are  rewritten  in  the  form: 


p  =  F,  -  u. 
*x  1  t 


D  =  F_  -  v, 

-  y  2  t 


so  that  differentiating  and  adding  we  find 


V  p  =  F,  +  F-  -  (u  +  v  ). 
r  lx  2y  x  y  t 


The  continuity  equation  (lc)  is  replaced  with  (5) .  In  the  global 
iteration  procedure  the  condition  (lc)  is  applied  in  obtaining 
appropriate  boundary  values  for  p  and  this  equation  is  satisfied 
through  the  convergence  of  the  pressure  solver. 

The  primary  differences  between  this  method  or  the  velocity- 
split  approach  for  subsonic  flows,  and  the  global  pressure 
relaxation  procedure  described  previously  for  thin  subsonic 
regions  are  two-fold:  (1)  the  pressure  gradient  in  the  momentum 
equations  are  treated  as  given  from  the  previous  iteration  and 
therefore  the  pressure  calculation  is  uncoupled  from  that  for  the 
velocities,  and  (2)  the  global  relaxation  includes  not  only  the 
pressure  but  the  velocities  as  well.  In  view  of  the  temporal 
gradients  appearing  in  (la),  (lb)  and  (5),  all  variables  are 
required  from  previous  iterations.  Solutions  cannot  be  obtained 
in  the  steady  state  mode  (At  =  °°)  discussed  in  the  previous 
section:  A  line  relaxation  procedure  coupling  (5)  with  (la) , 

(lb)  was  not  tested.  The  pressure  equations  were  solved 


independently  with  an  ADI  or  SOR  technique.  The  temporal 
gradients  in  (la)  and  (lb)  allow  for  smooth  passage  through 
separated  flow  regions. 

A  typical  solution  is  shown  in  figure  (5) .  These  PNS 
solutions  are  compared  with  full  Navier-Stokes  results  in 
reference  14.  The  agreement  is  good,  as  the  axial  diffusion 
effects  are  quite  small  even  in  the  regions  of  separated  flow. 
This  provides  another  example  where  steady  state  marching 
procedures,  with  global  pressure  relaxation,  can  be  used 
effectively  to  solve  the  PNS  equations  or  Navier-Stokes 
equations  where  the  influence  of  axial  diffusion  is  small.  A 
similar  approach  for  unbounded  subsonic  regions  is  being  tested 
with  the  velocity-split  method  mentioned  previously. 

Fig.  5 


5.  Example:  Cone  at  Incidence  -  Supersonic  Flow 

In  order  to  test  the  global  pressure  relaxation  procedure 

for  the  complete  compressible  PNS  equations,  the  supersonic 

flow  over  a  sharp  cone  at  incidence  has  been  considered. ^  The 

pressure  gradient  px  in  (la)  has  been  approximated  with  forward, 

central  and  backward  (single  sweep)  differences.  Lubard  and 
4  5  8 

Helliwell  and  Lin  and  Rubin  '  among  others  have  already 
investigated  the  latter  case.  Our  primary  interest  in  this 
regard  is  to  evaluate  the  limiting  value  of  (Ax)m^n  in  each  of 
the  cases  and  to  determine  whether  the  forward  difference 


16 


approximation  retains  the  effective  unconditional  stability 
predicted  in  the  previous  analysis. 

In  a  recent  study,  Schiff  and  Steger^  have  presented  a 

global  iteration  method,  but  not  with  the  intent  of  treating 

problems  where  upstream  interaction  or  axial  flow  separation 

are  important.  An  estimate  for  p^  is  obtained  by  an  initial 

backward  sweep  with  Ax  _>  (Ax)m^n«  Subsequently,  p^  is  treated 

as  known  from  the  previous  iterative  value  of  p.  In  addition, 

the  sublayer  approximation  presented  by  Rubin  and  Lin  (see 

reference  5)  is  aoplied  in  order  to  reduce  (Ax)  .  ;  i.e. ,  p 

min  x 

is  assumed  constant  across  the  subsonic  portion  of  the  boundary 
layer.  As  noted  earlier,  the  convergence  analysis  for  procedures 
in  which  px  is  treated  explicitly,  i.e.,  from  the  previous  sweep, 
would  indicate  that  these  global  relaxation  methods  are  unstable. 

In  reference  16  the  appearance  of  oscillations  is  noted  after 
four  global  iterations. 

A  second  important  feature  of  the  present  analysis,  that  is 
described  in  detail  in  reference  15,  is  the  application  of  the 
coupled  strongly  implicit  procedure  of  Rubin  and  Khosla^  for 
the  cross  plane  (normal  to  the  axial  flow  or  x  direction)  solution. 
This  is  considered  to  be  an  improvement  over  ADI  or  SOR  or  other 
splitting  techniques  as  there  is  an  immediate  coupling  of  all 
the  boundaries,  i.e.,  shock,  body,  lee  and  wind  planes.  The 
strongly  implicit  character  of  the  algorithm  also  appears  desirable 
for  capturing  imbedded  shock  waves  and  for  evaluating  secondary 
flow  separation  at  larger  angles  of  incidence.  In  addition, 


from  earlier  studies  convergence  rates  are  improved,  a  direct 

steady-state  solution  is  possible,  and  artificial  dissipation 

has  not  been  required  for  iterative  convergence. 

The  compressible  PNS  equations  are  given  by  an  expanded 

form  of  (1)  with  appropriate  energy  and  state  equations.  All 

12 

x  derivatives  are  backward  or  K-R  differenced  except  for  px 
which  is  forward  or  central  differenced  in  certain  cases.  A 
sublayer  approximation  is  not  assumed.  All  y  derivatives  are 
central  differenced.  The  trapezoidal  rule  is  used  for  the 
continuity  and  y-momentum  equations.  Boundary  layer-like  marching 
is  applied  in  the  x-direction;  the  normal  velocity  v  is  prescribed 
only  at  the  surface.  The  outer  value  of  v  is  coupled  with  and 
determined  by  the  Rankine-Hugoniot  conditions  at  the  shock. 

This  condition  also  provides  for  mass  conservation  in  the  shock 
layer. 

17 

The  strongly  implicit  algorithm  is  used  to  couple  the 
velocities  u,  v,  w  for  the  calculation  in  the  (y,£)  cross-plane. 
The  temperature  is  obtained  independently  from  a  similar  algorithm 
and  the  pressure  is  updated  from  the  normal  (y)  momentum  equation. 
The  strongly  implicit  algorithm  is  of  the  form: 


=  GMij 


Vi-lf j 


Vi, j+1 


where 


The  algorithm  (6)  is  scalarized  to  improve  computational  effi¬ 
ciency.  The  coefficients  GM..,  (■?..  0)  .  .  are  determined  from 

i j  i] 

recursive  formulas  whereby 


GM. 
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-  t<GtW 


GM 


(7) 
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etc.  The  boundary  conditions  on  all  surfaces  are  imposed 

either  in  (6)  or  (7) .  Symmetry  conditions  are  used  for  the 

wind  and  lee  planes  and  the  Rankine-Hugoniot  relations  are 

imposed  at  the  outer  shock  wave. 

Typical  solutions  are  given  in  figures  (6)  and  (7)  for  the 

shock  layer  thickness  and  the  heat  transfer  coefficient. 

-5  -+  -7 

Convergence  criteria  are  10  for  V  and  10  for  p  and  T.  The 

difference  grid  was  quite  coarse  in  the  cross  plane,  41x31  points. 

15 

With  finer  grids  the  accuracy  can  be  improved.  For  x/6  (x)  >_  2, 

where  6 (x)  is  the  shock  layer  thickness,  the  subsonic  layer 

thickness  YM  £  0-2  6.  Several  numerical  experiments  were  made 

in  order  to  estimate  the  values  of  (Ax)  .  .  The  results  are 

min 

given  in  Table  1.  Converged  results  were  obtained  with  central 
differencing  for  several  marching  steps;  however,  the  solutions 
were  oscillatory  and  extremely  inaccurate.  When  these  calculations 
were  continued  further,  the  instability  predicted  by  the 


Fig .  6 


Fig.  7 
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relation  (4)  soon  leads  to  divergence  of  the  solution.  With  the 
backward  and  forward  differences,  behavior  similar  to  that 
found  for  the  incompressible  equation  (1)  was  recovered  here. 

For  the  backward  px  calculations,  (Ax)m^n  appears  to  be  approaching 
an  asymptotic  limit  as  R  increases.  For  the  incompressible 

4 

problem,  at  R  -  10  ,  (Ax)  .  /y,.  ~  0.22,  as  opposed  to  the  value 

min  m 

0.25  obtained  for  the  cone  flow.  For  the  forward  p  calculations, 

x 

the  present  results  confirm  those  for  the  incompressible  equation; 
i.e.,  unconditional  stability.  Finally,  a  very  weak  stability 
condition  on  Ax  is  obtained  when  the  modified  "forward-difference" 
condition  is  applied,  see  Table  1.  As  seen  from  the  solutions 
presented  here,  this  approximation  is  quite  good  for  the  cone 
geometry.  Solutions  for  angles  of  incidence  up  to  45°  are 
presented  in  reference  15. 


R 

103 

Backward  p 

X 

0.2 

0.25 

Modified 
"Forward"  px 

10-3 

2  x  10"4 

Forward  px 

0 

0 

Central  px 

Unstable  for 
0(1) 

Unstable  for 
0(1) 

Table  1. 


Estimated  Values  of 
a  =  18°. 


(4xWyM!  M. 


7.95; 
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The  parabolized  Navier-Stokes  equations  have  been  considered 

for  subsonic  and  supersonic  flows.  It  has  been  shown  that  with 

single  sweep  marching  and  backward  differencing  for  all  axial 

derivatives,  including  the  pressure,  the  elliptic  influence 

is  numerically  suppressed  for  marching  steps  Ax  ^  (Ax)m^n  ~ 

where  y„  is  the  thickness  of  the  subsonic  zone.  For  subsonic 

boundary  layers  or  the  PNS  equations  with  supersonic  outer  flow 

conditions,  y  =  0(R  ^^)  and  therefore  (Ax)  .  =  0(R  . 

M  min 

For  triple  deck  regions,  y„  =  0(r”3//8)  or  (Ax)  .  =  0(R~3//8). 

M  mm 

For  problems  with  large  regions  of  subsonic  flow  (Ax)  .  =  0(1) 

mm 

and  large  truncation  errors  can  be  expected. 

If  global  relaxation  is  considered,  i.e.,  multiple  marching 
sweeps,  then  central  differencing  for  px  is  unstable,  but  forward 
Px  differencing  is  unconditionally  stable.  Forward  differencing 
also  has  the  desirable  property  of  allowing  for  a  complete 
pressure  coupling  with  the  velocities  and  a  free  surface  pressure 
interaction.  Therefore,  separation  regions  can  be  evaluated  with 
this  global  relaxation  method.  Examples  for  a  flat  plate  boundary 
layer  and  for  a  separation  bubble  have  been  presented. 

For  supersonic  outer  flows  the  global  relaxation  method 
requires  only  that  the  pressure  be  retained  from  the  previous 
iteration.  Therefore,  computer  storage  is  minimal.  With  large 
regions  of  subsonic  flow,  a  global  iteration  method  is  also 
presented;  however,  pressure  and  velocity  data  is  required  in 
this  procedure  and  therefore  computer  storage  will  be  increased. 


Solutions  are  presented  for  separated  channel  flow. 

Numerical  experiments  have  been  conducted  for  the  super¬ 
sonic  flow  over  a  cone  at  incidence.  A  coupled  strongly  implicit 
numerical  algorithm  was  applied  with  backward,  central  and 
forward  differencing  for  px.  The  solutions  confirm  the  analytic 
stability  results  for  the  incompressible  equations.  Backward 
differencing  is  conditionally  stable  {Ax  (Ax)m^n  2  ' 

central  differencing  is  unstable,  and  forward  differencing  is 
unconditionally  stable.  In  view  of  the  results  obtained  herein, 
we  conclude  that  for  flows,  with  thin  subsonic  layers,  forward 
differencing  for  the  px  term  leads  to  an  optimal  global  pressure 
relaxation  procedure,  with  free  pressure  interaction  and  minimum 
stability  limitations.  Global  relaxation  solutions  have  been 
obtained  for  subsonic  flows;  however,  optimization  of  such 
techniques  requires  further  study. 
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COMPARISON  OF  THE  LOCAL  SKIN  FRICTION  COEFFICIENT  FOR  A  FLAT  PLATE 
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FIG.  4  WALL  PRESSURE  VARIATION  THROUGH  THE  SEPARATION  ZONE 


Z-AXIS 


FIGURE  5.  STREAMLINE  CONTOURS  FOR  FLOW  IN  CHANNEL  :  R=103  (  FROM  REF. 
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FIGURE  7.  SYMMETRY  PLANE  HEAT  TRANSFER 
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for  thin  subsonic  boundary  layers  (yu  <<  1)  and  with  Ax  ■  0(yw),  stable  and 
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accurate  solutions  are  possible.  With  forward  differencing  of  the  axial 
pressure  gradient  the  procedure  can  be  made  unconditionally  stable;  a  global 
iteration  procedure,  requiring  only  the  storage  of  the  pressure  term,  has 
been  demonstrated  for  a  separated  flow  problem.  Solutions  for  incompressible 
boundary  layer-like  flows,  for  internal  flows,  and  for  supersonic  flow  over 
a  cone  at  incidence  with  a  coupled  strongly  implicit  procedure  are  presented. 
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